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Efficient  Numerical  Method  for  Computation 
of  Thermohydrodynamics  of  Laminar 
Lubricating  Films 


Prof.  Harold  G.  Elrod 
14  Cromwell  Court 
Old  Saybrook,  CT  06475 


Summary 


y  The  purpose  of  this  paper  is  to  describe  an  accurate,  yet 
economical,  method  for  computing  temperature  effects  in  laminar 
lubricating  films  in  two  dimensions.  The  procedure  presented  here  is 
a  sequel  to  one  presented  in  Leeds  in  1986  that  was  carried  out  for 
the  one-dimensional  case^J^ 

Because  of  the  marked  dependence  of  lubricant  viscosity  on  ^ 
temperature,  the  effect  of  viscosity  variation  both  across  and.  y 
along  a  lubricating  film  can  dwarf  other  deviations  from  MTdeal“— * 
constant-property  lubrication. 


"  C/#r  •"  " 

In  practice,  a  thermohydrodynamicsj  program  will  involve 


simultaneous  solution  of  the  film  lubrication  problem,  together 
with  heat  conduction  in  a  solid,  complex  structure.  The  extent  of 
computation  required  makes  economy  in  numerical  processing  of 
utmost  importance.  In  pursuit  of  such  economy,  we  here  use  techni¬ 
ques  similar  to  those  for  Gaussian  quadrature.  We  show  that,  for 
many  purposes,  the  use  of  just  two  properly  positioned  temperatures 
(Lobatto  points)  characterizes  well  the  transverse  temperature 
distribution . 
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1.  INTRODUCTION: 


The  purpose  of  this  paper  is  to  describe  an  accurate,  yet 
economical,  method  for  computing  temperature  effects  in  laminar 
lubricating  films.  The  procedure  presented  here  is  a  sequel  to  one 
presented  in  Leeds  in  1986. 1 

Because  of  the  marked  dependence  of  lubricant  viscosity  on 
temperature,  the  effect  of  viscosity  variation  both  across  and 
along  a  lubricating  film  can  dwarf  other  deviations  from  "ideal" 
constant-property  lubrication.  In  two  recent  papers  2-3  Khonsari  has 
summarized  the  growing  literature  concerned  with  the  present 
subject.  Consequently,  we  shall  not  undertake  a  survey  here,  but 
refer  only  to  those  articles  used  for  support  or  comparisons. 

In  practice,  a  thermohydrodynamics  A  program  will  involve 
simultaneous  solution  of  the  film  lubrication  problem,  together 
with  heat  conduction  in  a  solid,  complex  structure.  The  extent  of 
computation  required  makes  economy  in  numerical  processing  of 
utmost  importance.  In  pursuit  of  such  economy,  we  here  use  techni¬ 
ques  similar  to  those  for  Gaussian  quadrature.  We  show  that,  for 
many  purposes,  the  use  of  just  two  properly  positioned  temperatures 
characterizes  well  the  transverse  temperature  distribution. 

2.  NOMENCLATURE 

Single-underline  _  denotes  a  Legendre  coefficient. 
Double-underline  denotes  a  vector. 

Certain  matrices  are  defined  in  place  in  section  7. 

A  flow  vector,  defined  by  eq.  4.03 

B  flow  vector,  defined  by  eq.  4.04 

CP  specific  heat,  J/kg-C 

h  film  thickness,  m 

k  thermal  conductivity,  J/Csm 

l  subscript  for  lower  surface 

m  mass  flux  vector,  kg/smz 

p  pressure,  N/m2 

Pk  Legendre  polynomial  of  order  k 
t  time,  sec 

T  temperature,  deg  C 

u  velocity  of  fluid  in  x-direction,  m/s 
it  subscript  for  upper  surface 

v  velocity  of  fluid  in  y-direction,  m/s 
v  total  velocity  vector,  e>  u  +  eyv  +  e2w 

¥_  velocity  vector,  e x  u  +  e>v,  m7s  ~ 

w  velocity  of  fluid  in  z-airection,  m/s 
Wk  weight  for  ordinate  at  fk  ,  Lobatto  quadrature 
x  longitudinal  position  along  film,  m 
y  lateral  position  along  film,  m 

z  position  normal  to  film,  measured  from  midsurface 
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$  dissipation,  J/sm3  ,  defined  by  eq.  4.03 

p  viscosity,  Ns/m2 

£  fluidity  =  1/p,  m2 /Ns 

2z/h,  fractional  transverse  position 

p  density,  kg/m3 

3.  LOBATTO  INTERPOLATION  AND  QUADRATURE: 

In  general  terms,  the  problem  to  be  treated  is  one  of  three- 
dimensional  heat  convection.  A  numerical  solution  is  effected  by 
sampling  the  velocities,  pressure  and  temperature  over  a  chosen 
grid  of  points,  and  interlinking  their  values  by  algorithms  that 
incorporate  appropriate  physical  laws.  Generally  speaking,  the 
fewer  sample  points  required,  the  less  the  computational  effort. 
Consequently,  we  ask:  "What  are  the  best  cross-film  sampling 
positions?"  The  answer  is  provided  by  the  theory  of  orthogonal 
polynomials . 4  •  9 


Figure  1  shows  a  section  of  lubricating  film,  with  the  normal 
displacement,  Z,  scaled  to  be  -1  on  the  lower  wall,  and  +1  on  the 
upper.  At  certain  locations,  ^k  ,  we  intend  to  obtain  sample  values 
of  the  flow  variables  and  to  deduce  therefrom  certain  intermediate 
cross-film  values,  derivatives  and  integrals.  For  illustrative 
purposes,  consider  the  temperatures  Tk  =  T(fk).  These  are  assumed 
to  be  known  at  <f  =  l  and  <  =  -l,  and  at  N  intermediate  points. 

If  the  intermediate  sample  points  are  equispaced,  and  a 
polynomial  for  T(^)  is  passed  through  them,  the  excursions  of  such 
a  polynomial  between  points  can  become  unacceptably  large.  See,  for 
example.  Fig.  2,  where  a  tenth-degree  polynomial  is  used  to 
approximate  the  stepfunction  T(c<0)  =  1;  T(OO)  =  0  by  collocation 
at  eleven  equispaced  points.  On  the  other  hand,  if  a  tenth-degree 
polynomial  is  collocated  at  the  endpoints  (<r  =  -1;  <r  =  1)  and  at 
the  zeroes  (^k)  of  the  Jacobi  polynomial' a  >  Pb  1  •  1  (<)  ("Lobatto 
points"),  conformity  to  the  step  function  is  much  better.  Indeed, 
if  the  order  of  the  polynomial  is  increased,  equispace  interpola¬ 
tion  may  fail  to  converge  at  all,  whereas  the  latter  form  of 
interpolation  becomes  progressively  better. 


Not  only  do  the  Lobatto  points  serve  well  for  interpolation 
by  high-order  polynomials,  but  they  also  serve  for  accurate 
numerical  integration3  .  It  can  be  shown  that  N  such  internally- 
selected  points  permit  exact  numerical  integration  of  a  polynomial 
of  order  2N+1  over  the  range  -l<f<l.  Thus: 


—  -  A 

<*>  For  N  interior  points,  the  Lobatto  locations  (<\.)  are  at  the  ■  —  . . 

zeroes  of  dP*»  i  (  i  )  /d<  =  N  (N+l )  P*- 1  *  •  1  ( <T  )  ,  where  Pn  {<  )  is  the  Legendre  _ 
polynomial  of  order  N. 
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Availability  Ooda# 


N 


Diet 


avail  and/or 
Spaolal 


□O 


1 

[2.01]  JT(f)df  =  EwuTk 
-1 

The  Lobatto  locations,  £k ,  and  weight  factors,  Wk ,  are  to  be  found 
in  Abramowitz  and  Stegujfc .  For  N=2,  we  have: 

Location 
-1 

-1/J5 
1/J5 
1 

Note  that,  in  contrast  to  Gaussian  quadrature,  Lobatto' s  technique 
incorporates  endpoint  values  which,  in  our  case,  are  known,  and 
might  as  well  be  used. 

The  approach  taken  in  this  paper  is  to  form  partial  differen¬ 
tial  equations  in  x,y  for  the  Lobatto-point  temperatures,  with  the 
transverse  temperature  distribution  given  by  a  collocated  polynomi¬ 
al.  The  "fluidity",  £  =  1/p  is  also  collocated  to  its  Lobatto-point 
values ,  Ck  =  i (Tk ) . 

Let  us  turn  now  to  the  analytical  development  which  incor¬ 
porates  these  ideas. 

3.  BASIC  THERMOHYDRODYNAMIC  EQUATIONS: 

The  equations  to  be  used  are  for  laminar  lubrication  with  a 
fluid  possessing  constant  density  and  constant  thermal  properties. 
The  momentum  equation  is: 

[4.01]  (a/dz) (pOv/az)  =  7p 

the  pressure  being  independent  of  "z" . 

And  the  energy  equation  is: 

[4.02]  pCP  (DT/Dt)  =  k(d*T/az2)  +  4  with: 

[4.03]  $  =  p(av/az)2 

Here  the  viscosity,  p,  may  depend  markedly  on  the  local  temperature 
T.  In  addition,  the  conservation  of  mass  (or  volume)  must  be 
satisfied;  i .  e. , 

[4.04]  7-v  =  0 

Several  investigators7’8  have  justified  the  use  of  these 
equations  by  showing  that  the  effects  of  variations  in  density, 
specific  heat,  etc.,  of  the  lubricant  are  quite  secondary  to  those 
resulting  from  changes  in  its  viscosity,  which  usual) v  varies 
substantially  with  temperature. 

h 


The  boundary  conditions  to  be  satisfied  by  eqs.  [4.01  -  4.03] 
are  as  follows.  The  liquid  velocities  must  conform  to  the  upper 
and  lower  surface  motions,  and  the  liquid  temperature  to  the 
surface  temperatures.  At  the  film  peripheries,  incoming  liquid 
temperatures  correspond  to  ambient  conditions,  whereas  outgoing 
temperatures  derive  from  the  interior  where  they  originate. 

5.  VELOCITIES.  REYNOLDS  EQUATION: 

To  compute  the  velocities,  it  is,  for  several  reasons,  conven¬ 
ient  to  use  the  fluidity  expressed  as  a  series  in  Legendre 
polynomials.  Thus: 


[5.oi]  c  =  SLk Pk  (<r) 


As  previously  stated,  the  coefficients  £>•  are  here  determined  by 
collocation  of  the  fluidity  at  the  Lobatto  points,  where  its  values 
are  determined  by  the  Tu  .  A  double  integration  of  eq.  [5.01]  then 
yields  for  V  =  ex  u  +  e>  v . 

<  <T 


[5.02] 
where : 
[5.03] 
and : 


V  =  Vi  +  A  Si  +  B  Sit  d{ 

-1  -1 

1  1 

A  =  [Vo  -  Vi  -  B  m<l<r]/[.rcdn 

-1  -1 


[5.04]  B  =  (h/2) 2  Vp 


To  obtain  the  lineal  mass  flux,  eq.  [5.02]  must  be  integrated 
again.  The  result  is: 


[5.05]  m/p  =  (Vu  +  V  l  )  ( h  /  2 )  -  (h/3)£_iA  -(2/3)B(£_o  +  2^/5) 

This  result  is  independent  of  the  number  of  terms  in  the  series 
[5.01].  When  it  is  inserted  into  the  mass  continuity  equation: 

[5.06]  dh/at  +  7* (m/p)  =  0 

the  following  generalized  Reynolds  equation  results: 


[5.07]  <7*  t  P  h3  9p  =  6  ( Vu  +  Vi  )  *7h  +  12  (ah/at)  -  27*  (Cl  /Co  )h(Vo  -  V.  ) 

where : 


[5.08]  Cp  =  L®  +  0.4Lz  “  (La  )2  /(3£j>  ) 


In  form,  eq.  [5.07]  differs  from  the  standard  Reynolds  equation'* 
only  through  its  last  term. 
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Now  we  wish  to  express  our  differential  equations  using  the 
fractional  gap  position  as  transverse  coordinate  across  the  film. 
Thus : 

[5.09]  Z  =  2z/h 

where  "z"  is  measured  from  the  midsurface.  Taking  "x"  as  a  typical 
lateral  coordinate  and  "u"  as  a  typical  variable,  we  note  that: 


[5.10]  (du/ax)y,z  =  (au/ay)y,£  -  ^  ( aiog (h) /ax) y  (au/a«r ) x , y 

The  transverse  velocity,  w,  can  be  found  from  mass  continuity 

via : 

[5.11]  aw/az  =  -(au/ax  +  av/ay)  = 

But: 


[5.12]  (r*V)z  =  -  ^log(h)  •  (dV/dZ) 

Substituting  [5.12]  into  [5.11]  and  integrating,  we  get: 

z 

[5.13]  W  =  wL  -  Xl9<r*v  -  ^9iog(h)  •  (av/a-r)  Idz 

-h/2 

Integration  by  parts  gives: 


Z 

[5.14]  w  =  WL  +  Vi  *?(h/2)  +  £V*<?(h/2)  -  VZ  *  [  (h/2 )  jVd<r  ] 

-1  ~ 


6.  TEMPERATURE  EQUATION: 

To  obtain  a  convenient  differential  equation  for  the  tempera¬ 
ture,  we  rewrite  eq.  [4.02]  as: 

[6.01] 

aT/at  +  u(aT/ax)y.^  +  v(aT/ay)x,f  +  (2/h)  (aT/a<r)x.v  (w  -<rv*^(h/2)  i  = 

( 4tc/h2  )  (  a2  T/dZ2  )  x  ,  y  +  4/(pCp) 

To  treat  the  cross-velocity  term,  we  note  that  kinematics  gives: 
[6.02]  WL  +  V*9(h/2)  =  -a(h/2)/at 

Substituting  [5.14]  and  [6.02]  into  [6.01],  we  get: 


<r 

[6.03]  aT/at  +  V7T  -  (l/h)  (aT/d£)  [  (i+<r)  (ah/at)  +  ?‘hjvd<:]  = 

-l 

(4K/h2 )  (a2T/a<r2 )  +  4/ (pcP ) 
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All  partial  derivatives  in  the  foregoing  equation  are  in  x,y,f,t 
space.  There  is  one  such  partial  differential  equation  16.03]  in 
x , y , t  for  each  Tk . 


7.  NUMERICAL  PROCEDURES: 

It  is  convenient  to  treat  the  Tk  (xi , yj  ) ,  tk  (xt  , yj )  as  vector 
matrices,  and  to  use  matrix  notation.  We  have: 

[7.01]  T  U)  =  vt„  P„  U) 

[7.02]  t  U)  =  ZL«  P«  U) 

Therefore : 

[7.02]  Tk  =  ZT„  Pa(<Tk)  or:  Tk  =  ECk  n  Tk  or:  T  =  CT  and:  T  =  C~lT 
Also,  then: 

[7.03]  t  =  CL  and:  L  =  C*'  t 

Successive  differentiation  of  [7.01]  enables  us  to  write: 

[7.03]  aT/a<T  =  DT  and:  a2T/d^2  =  ET 

Successive  integration  of  [7.02]  permits  us  to  write 

[7.04]  v  =  VL  +  AFC  +  BGt  and: 

'C 

[7.05]  /Vd^  =  Vl ( 1+<T )  +  ARC  +  BSC 
-1 

To  preserve  numerical  stability,  backward  differences  are  used 
for  the  convective  terms.  Thus,  at  the  point  (i,j)  we  write: 

[7.06]  u(aT/dx)  =  abs  (ui ,  j  )  [Ti ,  j  -  Ti-»x,j];  gx  =  sgn(ui.j) 

with  a  similar  expression  for  v(aT/ay). 

Now  form  the  diagonal  matrix,  Ed,  from  the  diagonal  terms  of 
E,  and  let  E0  =  E  -  Ed  .  With  these  definitions,  the  matrix  finite- 
difference  equation  corresponding  to  [6.03]  is: 

[7.07]  [1/At  +  abs (u) / Ax  +  abs(v)/Ay  -  ( 4k /h2 ) Ed  ]  T"p w  i  .  ,  = 

Tt.j/At  +  t  abs  (u)  Ti-ox  .  j  I /Ax  +  (  abs  (v)  Ti  , j-9y I  /Ay  +  (4ic/h2)Eo  Ti  ,  i  + 

( 1/h)  (V*hjVd<r)D  Ti  .  j  +  4/pCp 
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Here  Tnew  denotes  T(t+At),  u,  v,  $  and  T  belong  to  the  same  trans¬ 
verse  row  (<Tk),  and  the  integral  of  velocity  is  from  the  lower  wall 
to  fk  . 


The  explicit  form  of  the  difference  equation  [7.07]  is  stable 
for  a  large  At,  such  as  was  used  to  generate  the  steady-state 
results  presented  in  this  report. 

8.  COMPUTATIONAL  RESULTS: 

Validation  of  the  computer  program  was  attempted  in  various 
ways.  Comparisons  were  made  with  the  one-dimensional  calculations 
of  Hunter  and  Zienkiewicz10  ,  Dowson  and  Hudson7  and  of  Hahn  and 
Kettleborough8 ,  who  all  used  the  physical  properties  tabulated 
below . 

p  =  0.13885  expl-p  (T-Ta.bient  )  I  ,  Ns/m*  ;  [3  =  0.045 
pCP  =  1.7577  J/m:l  C 
k  =  7 . 306  E-08  m* /s 

Tambtent  =  0  (used  as  reference,  only) 

These  same  properties  are  used  for  all  examples  of  this  report, 
except  that  [3  =  0,  instead  of  0.045,  for  the  constant-property 
calculations . 

Dowson  and  Hudson  chose  the  following  bearing  characteristics. 
L  =  0.18288  m;  W  infinite 
Uu  «  31.96  m/s;  UL  =  0 

hi  =  1.8288  E-04  m;  h2  =  0.9144  E-04  m 

Figures  3  and  4  compare  our  results  with  theirs<b>.  Agreement  is 
certainly  satisfactory,  though  not  perfect.  Dowson  and  Hudson  used 
Ax  =  L/20  and  a  constant  Ay  =  h2/20,  whereas  we  employed  Ax  =  L/30 
and  N  =  8.  In  these  comparisons,  as  in  others  made,  the  reasons  for 
the  observed  small  discrepancies  are  difficult  to  assign.  There¬ 
fore,  it  was  decided  to  write  the  present  program  for  arbitrary  N, 
and  to  use  N  =  8  as  a  standard  against  which  to  test  more  approxi¬ 
mate,  but  faster  versions  employing  N  *  2  or  3.  To  assist  others 


<b)  All  dotted  curves,  except  those  in  Fig.  2,  were  collocated  at 
the  Lobatto  points,  and  interpolated  with  an  auxiliary  plotting 
program;  i.e.  the  dotted  curves  are  not  everywhere  in  complete 
agreement  with  the  corresponding  Legendre  polynomial  expansions. 
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who  may  wish  to  make  comparisons  with  their  own  programs,  we  list 
some  results  for  the  above  problem  in  the  Appendix. 

In  Fig.  3,  the  great  reduction  in  load  capacity  due  to 
lubricant  warmup  is  manifest.  Of  course,  in  engineering  design 
practice  some  allowance  for  this  effect  is  made  by  assigning  some 
lowered  constant  viscosity  corresponding  to  an  estimated  tempera¬ 
ture  rise. 

Figure  4  shows  the  convergence  with  N  of  results  obtained 
using  the  present  program.  There  is  no  perceptible  difference 
between  the  results  for  N  =  5  and  those  for  N  =  8.  Somewhat 
fortuitously,  the  results  for  N  =  2  also  almost  coincide  with  the 
"true"  curve.  Figure  5  shows  the  convergence  of  the  temperature 
distributions.  The  Lobatto-point  results  for  N  =  2  and  N  =  3  are 
compared  with  the  curves  for  N  =  8  at  the  bearing  exit,  and  at  the 
halfway  point. 

Reverse  convection  has  been  a  source  of  difficulty  for  some 
investigators.  Accordingly,  we  show  in  Figs.  6  and  7  some  calcula¬ 
tions  for  a  high  f  ilm- thickness  ratio  of  4.  Note  the  rapid 

variation  of  temperature  in  the  bearing  inlet  -  incoming 

temperatures  are  taken  at  an  entrance  value  of  0,  whereas  tempera¬ 
tures  in  the  backflow  region  originate  from  a  region  where  viscous 
heating  has  occurred.  As  shown  in  Fig.  1,  Lobatto  interpolation  is 
particularly  suited  to  cope  with  such  variation.  The  two-tempera¬ 
ture  version  (N  =  2)  of  our  program  performs  surprisingly  well. 

All  of  the  foregoing  calculations  are  for  a  one-dimensional 
bearing,  run  with  the  program,  however,  as  a  two-dimensional  wide 
bearing.  The  same  cases  were  also  solved  earlier  by  a  one-dimen¬ 
sional  procedure*  embodying  Lobatto-point  methods.  But  when 
extended  to  two-dimensional  problems,  that  earlier  procedure  proved 
to  be  persistently  unstable.  To  deal  successfully  with  the  added 
dimension,  the  present  method  was  devised.  Figures  8,  9  and  10  show 
some  sample  "true"  two-dimensional  results  obtained  with  a  square 
slider  bearing. 

We  intend  to  continue  this  investigation  by  engaging  in  some 
parametric  studies  and  by  coupling  the  new  technique  to: 

inlet-groove  temperature  distribution 
cavitation 

heat  conduction  in  the  bearing  body 
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APPENDIX: 


DISTRIBUTION 

FOR  DOWSON-HUDSON  PROBLEM 

x/L 

pressure,  N/m2 

0 

0 

.  1 

3545479 

.2 

6434485 

.  3 

8755968 

.  4 

1 . 051984E+07 

.  5 

1 . 167  347E+07 

.6 

1 . 2097 84E+07 

.7 

1 . 158809E+07 

.  8 

9815331 

.9 

6259130 

1 

0 

TEMPERATURES  OBTAINED  FOR  DOWSON-HUDSON  PROBLEM 

X/L 


> 

0_ 

0.2 

0 . 4 

0 . 6 

0 . 8 

1.0 

-1 
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Figure  6.  Temperature  Profiles  at  Exit  and  Halfway 
Infinite-Width  Slider,  Variable  Viscosity,  h» /hi  =  2 
N=8  {...);  N=3  (*) ;  N=2  (X) 
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Figure  7.  Pressure  Profile  for  Infinite-Width  Slider 
Variable  Viscosity,  hi /ha  =  4 
N=8  (...);  N=2  (X) 
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Variable  Viscosity,  h» /hi  =  2 
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